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Abstract 

We report on the retrieval of directional sea-roughness (the full directional mean square slope, including 
MSS, direction and isotropy) through inversion of Global Navigation Satellite System Reflections (GNSS- 
R) and SOlar REflectance Speculometry (SORES) data collected during an experimental flight at 1000 m. 
The emphasis is on the utilization of the entire Delay-Doppler Map (for GNSS-R) or Tilt Azimuth Map (for 
SORES) in order to infer these directional parameters. Obtained estimations are analyzed and compared to 
Jason- 1 measurements and the ECMWF numerical weather model. 

Keywords: GNSS-R, GPS, Galileo, Speculometry, SORES, sea roughness, Directional Mean Square 
Slope, Delay-Doppler Map. 



1. Introduction 

The use of Global Navigation Satellite System 
(GNSS) signals reflected by the sea surface as a 
remote sensing tool has generated considerable at- 
tention for over a decade 1121 1131 . Among sev- 
eral applications, two classes have rapidly emerged 
in the community: sea-surface altimetry, which 
aims at retrieving the mean sea level like classical 
radar altimeters do, and sea-surface scatterometry or 
"speculometry" (see below for a justification of this 
neologism) for the determination of sea roughness 
and near surface wind. This paper addresses the lat- 
ter application. 

Inferring sea roughness from GNSS-R data re- 
quires (i) a parametric description of the sea surface, 
(ii) an electromagnetic model for sea-surface scat- 
tering at L-band and (iii) the choice of a GNSS-R 
data product to be inverted. In the literature, there 
is quite an agreement on the two first aspects. It 
has been recognized that the scattering of GNSS 
signals can be modeled as a Geometric Optics pro- 
cess (GO), where the fundamental physical process 
is the scattering from mirror-like surface elements. 
This is the reason why we use the term "speculome- 
try", from the latin word for mirror, speculo. There- 
fore, the most important feature of the sea surface 



is the statistics of facet slopes at about the same 
scale as the electromagnetic wavelength (A). This 
is described by the bi-dimensional slope probabil- 
ity density function (PDF). Under a Gaussian as- 
sumption, three parameters suffice to fully define the 
sea-surface PDF: the directional mean square slope 
DMSSa, which results from the integration of the 
ocean energy spectrum at wavelengths larger than A. 
The symbol DMSSa englobes the three parameters 
defining the ellipsoidal shape of the slope PDF: scale 
(total MSS), direction (Slope PDF azimuth) and 
isotropy (Slope PDF isotropy). The GNSS-R scat- 
tering model proposed by Zavorotny and Voronovich 
in 1171 is based on GO, and is, to date, the refer- 
ence model for the GNSS-R community. While for 
the purposes of specular scattering the sea-surface 
roughness is parametrized by the directional mean 
square slope in a direct manner, DMSSa is rarely 
emphasized as the geophysical parameter of interest. 
Instead, most authors prefer to link sea roughness 
to the near surface wind vector, which is thought to 
be more useful for oceanographic and meteorologi- 
cal users, but misleading. Indeed, this link requires 
an additional modeling layer and is an extra source 
of error. For instance, a wind-driven sea spectra is 
not suitable for inferring sea surface DMSSa when 
swell is present or the sea not fully developed. The 
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connection between DMSSa and wind is thus mod- 
ulated by other factors (e.g., swell, fetch and degree 
of maturity). 

Usually, for technical reasons, the product in- 
verted in GNSS-R speculometry is a simple De- 
lay Waveform, that is, a ID delay map of the re- 
flected signal amplitude. Using a single GNSS 
emitter, the wind speed can be inferred assuming 
an isotropic slope PDF (i.e., the PDF's shape is a 
circle) l8l l3l fTTl . Attempts have also been made 
to measure the wind direction by fixing the PDF 
isotropy to some theoretical value (around 0.7) and 
using at least two satellites reflections with differ- 
ent azimuths 1181 l9l. As investigated in the frame 
of the ESA OPPSCAT project (see |lj and H51 ). 
it is nonetheless possible to work on a product of 
higher information content: the Delay-Doppler Map 
(DDM), a 2D delay-Doppler map of the reflected 
signal amplitude. The provision of an extra dimen- 
sion opens the possibility to perform the full estima- 
tion of the DMSSa- In Q, Elfouhaily et al. devel- 
oped a rapid but sub-optimal method based on the 
moments of the DDM to estimate the full DMSSa: 
this approach neglects the impact of the bistatic 
Woodward Ambiguity Function modulation of the 
Delay Doppler return. 

The present paper was motivated by a recent ex- 
periment conducted by Starlab for the demonstration 
of GNSS-R altimetry. The altimetric aspects are re- 
ported elsewhere 1151 . We note that the configura- 
tion of the flight was not optimized for speculom- 
etry: from 1000 m altitude, the sea-surface reflec- 
tive area is essentially limited by the PRN C/A code, 
and the glistening zone is coarsely delay-Doppler 
mapped. In addition to the GNSS-R experiment, 
high resolution optical photographs of sun glitter 
were also taken, providing the SORES dataset (SO- 
lar REflectance Speculometer). Since the classic 
study of Cox and Munk |4|, it is well known that 
sea-surface DMSSopt can be inferred from such 
data. The availability of optical photographs thus 
provided us with an extra source of colocated infor- 
mation. Because there is a strong similarity between 
products —DDM for GNSS-R and the Tilt Azimuth 
Map (TAM) for SORES— and models —GO— the 
same inversion methodology can be applied to both 
datasets. 

The goal of the paper is to investigate the full 
exploitation of the bidimensional GNSS-R DDM 
and SORES TAM products to infer the set of three 
DMSSa parameters. The driver of the study has 
been the exhaustive exploitation of the information 



contained in those 2D products. Consequently, the 
proposed approach relies on a least-square fit of 
speculometric models to datasets. We first describe 
in details the collected datasets and the associated 
pre-processing. Then, we present the speculometric 
models needed to infer data together with the inver- 
sion scheme. Finally, we provide the estimation re- 
sults and discuss their coherence with other sources 
of data. 



2. Dataset collection and pre- 
processing 

The campaign took place Friday September 27th 
2002 around 10:00 UTC, along the Catalan coast 
from Barcelona (Spain) up to the French border. An 
aircraft overflew the Jason- 1 track 187 at 1000 m 
along 150 km and gathered 1.5 hours of GPS-R raw 
signals (see |16| for more details). Since it would 
have been computationaly too expensive to process 
the full track, it was divided into 46 10-second arcs 
(each spanning roughly 500 meters), sampled every 
50 s (see Figure Q). The first arc started at GPS 
Second Of the Week 468008.63. The aircraft kine- 
matic parameters were kept close to the nominal val- 
ues specified in the mission plan: altitude=1000 m, 
speed=50 m/s and heading from North=30°. We 
have selected three GPS satellites in optimal view 
during the experiment whose elevation and azimuth 
are given in Figure|2] 

The raw GNSS-R data were recorded using the 
GPS reflection equipment provided by ESA. Specif- 
ically, the GPS direct and reflected signals were 1- 
bit sampled and stored a rate of 20.456 Mbit/s. The 
pre-processing step consisted in performing a delay- 
Doppler Pseudo Random Noise (PRN) code de- 
spreading to coherently detect the direct signal (from 
GPS emitter) and the reflected signal (scattered by 
sea-surface). We used the Starlab in-house soft- 
ware to produce three DDMs time-series (one per 
PRN), sampled into 46 arcs of 10 seconds each. The 
general processing strategy was to track the delay- 
Doppler of direct signal and then compute DDMs for 
both direct and reflected signals. These DDMs ac- 
tually represent the filtered electromagnetic field of 
incoming signals, as processed with delay-Doppler 
value slightly different from those corresponding to 
the specular point. The coherent integration time 
was set to 20 ms to ensure a Doppler resolution of 
50 Hz. The delay map spanned 80 correlation lags 
(i.e. +/- 1.95 /is) with a lag step of 48.9 ns, while the 
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Figure 1: Map of the aircraft track divided into 46 10-second arcs. 



Seconds after SOW=4 
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Seconds after SOW =4 



Figure 2: Elevation and Azimuth of three GPS satellites in view during the 46 10-second arcs. 



Doppler range spans -200 Hz to 200 Hz with a step 
of 20 Hz. Incoherent averaging was applied to each 
arc (the accumulation time was set to 10 s). This 
process aimed at reducing both thermal and speckle 
noise by a factor of a/500. At the end, the GNSS- 
R product for one PRN and one arc was an average 
amplitude delay-Doppler map of size 81x21. 

The SORES photographs were taken from time 
to time along the track when the roll, pitch and yaw 
angles of the plane were negligible. The camera 
was a Leica dedicated to aerial photography. An 
inertial system (by Applanix) provided the time- 
tagged position for each snapshot. The film was a 
panchromatic Aviphot Pan 80. The focal length was 
15.2 cm, and the photographic plate had an area of 
23x23 cm 2 . The aperture angle was consequently 



74.2°. The observed area was a square with area 
1. 1 24 x 1.124 km 2 . The exposure time was fixed to 
1/380 and aperture to f/4. The silver photographs 
were scanned and digital images were averaged to 
400x400 squared pixels, in order to be easily pro- 
cessed. 

3. Models and inversion 

3.1. Directional mean square slope 

As discussed above, forward scattered signals at 
short wavelengths (optical but also L-band) are 
mostly sensitive to the specular scattering mecha- 
nism. Therefore, the strongest sea-surface signa- 
ture in received signal is expected to be due to facet 
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slope statistics at the relevant scales. The 2D slope 
probability density function (PDF) Gaussian model 
is given by 



V{s) 



(1) 



where s' = [dQ/dx d(/dy] stands for the vector 
of directional slopes in some frame xy and K is 
the matrix of slope second order moments. The xy 
frame mapped on sea-surface is defined as follows: 
it is centred on the specular point and has its x axis 
aligned with the Transmitter-Receiver line. Mean- 
square slopes along major and minor principal axes 
are often referred to as MSS up-wind (mss u ) and 
MSS cross-wind (mss c ) respectively. The K matrix 
is then obtained via a simple rotation, 



K = Em 



mss u 
mss c 



■ R- 



(2) 



where R^ is the usual rotation matrix of angle ip, 
the angle between the x-axis and the slope principal 
axis. 

Thus, mss u , mss c and ip are the three geophys- 
ical parameters we wish to estimate. They can be 
thought of as the three parameters of an ellipse (see 
Figure [3} representing the slope PDF mapped on 
sea-surface. In the following, we will consider the 
equivalent set of parameters: 

• Total MSS, defined as: 2y/mss u .mss c - This 
magnitude is actually proportional to ellipse 
area and can be interpreted in terms of wind 
speed, based on a particular wind-driven sea- 
surface spectrum like the Elfouhaily's spec- 
trum 1 6 1 . 

• Slope PDF azimuth (SPA), defined as the 
direction of semi-major axis with respect 
toNorth. As shown by Figure [3] this angle is 
7r + $ — ip, if $ i s the satellite azimuth from 
North. 

• Slope PDF isotropy (SPI), defined as 
mss c /mss M . When SPI=1, the slope PDF 
is isotropic and the glistening zone is circular. 
Low values of SPI indicate a highly directive 
PDF. Typically, SPI is expected to be around 
0.65 for well developed sea-surface. 

3.2. GNSS-R speculometric model 

The classical GNSS-R bistatic radar equation 1171 
links the average GNSS-R power return to sea- 
surface slope PDF. Processing the raw signal with 



various delay-Doppler values (r, /), a DDM is com- 
puted whose theoretical expression is: 



dxdy ■ 



Cj r 



X 2 ( T m -T c -T, f m - f c - f) (3) 

where G r is the receiver antenna pattern, R t and R r 
are the distances from generic point on sea-surface 
to transmitter and receiver, <r° is the reflectivity, \ is 
the Woodward Ambiguity Function (WAF, see 1131 ). 
T m (i, y) and f m (x, y) are the delay-Doppler map- 
ping on sea-surface and t c and f c are delay-Doppler 
centers. To first order, the reflectivity is proportional 
to the slope PDF: 



-q.r 



-jty 



(4) 



where (q x , q y , q z ) is the scattering vector and 1 1Z | 2 = 
0.65 is the specular Fresnel coefficient. 

The presence of thermal noise biases the value of 
average power return. Hence, the average amplitude 
of the DDM can by modeled by 



A(r, f) - y/a.P(T,f) + b, 



(5) 



where b stands for the bias in power. In particular, 
this effect is visible in the early-delay domain of the 
DDM: for delays lower than one-chip, the DDM am- 
plitude has a no null value, often called "grass level". 
As we do not have a calibrated model an overall scal- 
ing parameter a is also needed in the model. 

To sum up, the model features three parameters 
of interest and four "nuisance parameters": 

• the DMSSa, characterizing the Gaussian 
slope PDF: total MSS, isotropy (SPI) and az- 
imuth (SPA), 

• the DDM delay-Doppler centers: t c and f c , 

• overall scaling parameter: a, 

• grass level: b. 

Other parameters required to run the forward model 
are recalled in Tabled 
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Figure 3: Sketch of the slope PDF and related frames. 



Geometry 


Aircraft: Altitude, speed and heading provided at 1 Hz 
Satellite: Elevation and azimuth provided at 1 Hz 


Instrument 


Antenna Pattern: measured in anechoic chamber 
Band: GPS LI (19 cm) 
GPS Code: C/A 


Processing 


Integration Time: 20 ms 
Accumulation Time: 10 s 

Doppler span: [ -200 Hz , 200 Hz ], 20 Hz step 

Delay span: [-40 samples, 40 samples], 1 sample step 



Table 1 : Overview of the parameters necessary for running the DDM forward model. 



3.3. SORES speculometric model 

To date, results derived from the glitter-pattern of re- 
flected sunlight as photographed by Cox and Munk 
in 1951 remain the most reliable direct measure- 
ments of wind-dependent slope statistics. As ex- 
plained in their well-documented report |5 1, the sea- 
surface can be gridded with a Tilt (/3) Azimuth (a) 
Mapping of the small facet slopes. These are just a 
polar parametrization of the (s x ,s y ) slopes discussed 
in the previous section: 



cos a 
sin a 



tan/3, 
tan/3. 



(6) 



The link between the small facet slope PDF V and 
the intensity in the photograh I m is given by 



I m (a,0)=A o -V(a,(3)- f(a,(3, 
+K-I b (a,P), 



(7) 



where Ib(a,(3) is the intensity of the picture back- 
ground (i.e. far from the sun glint), K and Aq are 



multiplicative constants and / is a transfer function, 

ft r a\ P(l - p) 3 sm0cos 3 /x 

f(a,p,(p) = 3-5 , (8) 

cos^ p cos Ul 

with <fi: sun elevation, p: coefficient of reflection and 
/1, lu: two angles shown on Figure 

The pixel intensity on the image comes princi- 
pally from the additive contribution of sunglint and 
reflected skylight. The sunlight scattered by parti- 
cles beneath the sea surface is assumed negligible 
and is not considered here. A model has been de- 
veloped to remove reflections of sky radiance from 
the glint. The approach consists in considering each 
sea surface facet specular because, for a given loca- 
tion of the receiver, there always exists a "source" 
in the sky satisfying the specular reflection condi- 
tion. Let's consider the cell (on, f3i) of the TAM. It 
corresponds to the slope components required to re- 
flect the solar rays onto the camera. The radiance of 
the sea surface due to reflected skylight in the cell 
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Figure 4: Geometry of the SORES experiment. 



Photograph 41 




meters 

Figure 5: Tilt Azimuth Mapping overlaid on a SORES photograph. 

(a i; Pi) can thus be simply modeled by the integra- 3.4. Inversion Scheme 
tion of intensity over all the azimuths a and tilts (3 

except the azimuth cti and the tilt # of the corre- Inversion was performed through a minimization 
sponding cell. of the root mean square difference between model 

and data products (i.e., DDMs for GNSS-R and 
TAMs for SORES). Numerical optimization was 
carried out with a steepest-slope-descent algorithm 
(Levenberg-Marquardt type adjustment). The three 
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DMSSa as well as nuisance parameters (r c , f c and 
a for the DDMs and A Q and K for the TAMs) 
were jointly estimated in an iterative manner: nui- 
sance parameters (as a first step) and DMSSa (as a 
second step) were successively optimized, repeating 
this two-step sequence until convergence. Figure [6] 
gives qualitative examples of fit results for DDM and 
TAM. 

4. Results and discussion 

Figure shows total MSS, Slope PDF Azimuth 
and Slope PDF Isotropy along the aircraft track be- 
tween latitudes 41.2° and 42.2°, as estimated by 
SORES and GNSS-R. Other sources of information 
are also shown for comparison: 

• total MSS for C- and Ku- bands, derived from 
the Jason- 1 o° measurements at 1 Hz sam- 
pling (7 km) and 20 km resolution, 

• wind direction provided by the ECMWF nu- 
merical weather model, with accuracy of 20°, 
and 

• swell direction derived from a spectral analy- 
sis of SORES images (the periodic pattern of 
long waves is indeed clearly observed on the 
photographs). 

4.1. Total MSS 

The total MSS has been plotted in log-scale in order 
to compare different frequency measurements more 
easily. The common trend for all bands is the in- 
crease of slope variance with latitude until a relative 
plateau is reached. Measurements of PRN08 and 
24 show good agreement while PRN10 seems to be 
somewhat up-shifted. As expected, we observe that 
the level and dynamic of MSS decrease with longer 
wavelength: Optical, Ku, C and L band, in this or- 
der. Nevertheless, the level and dynamic range of 
GNSS-R plots (especially PRN10) seem a bit large 
for L-band measurements, when compared to C- 
band. Note however that Jason- l's MSS have been 
obtained through the relationship MSS — n/a°, 
K being an empirical parameter accounting for cal- 
ibration offsets. Unfortunately, the uncertainty on 
k makes the absolute levels of Jason- l's plots very 
doubtful. Here, as an illustration purpose, we have 

* The Jason- 1 wind speed is derived from Ku-band <r° and signil 
accuracy is 2 mis. 



set k=0.45 and k=0.95 for C- and Ku-band respec- 
tively. 

4.2. Slope PDF Azimuth 

Using a single DDM, the estimation of SPA is de- 
generate in two particular cases: when the transmit- 
ter is at zenith or when the receiver moves torwards 
the transmitter |2|. In these two cases, the Delay- 
Doppler lines that map the glistening zone are sym- 
metric around the receiver direction. Hence, one 
cannot distinguish between a slope PDF and its mir- 
ror image about the receiver direction. Here, PRN08 
has its elevation comprised between 74 and 83 de- 
grees. It is then very likely that the SPA estimated 
for this PRN is degenerate. For this reason, we have 
added on the plot the mirror image of the SPA about 
the receiver direction (30°). We also note that the az- 
imuth of PRN 10 decreases down to 230° at the end 
of the track, quite close from 210°, the complemen- 
tary of the receiver's direction. 

According to ECMWF data and SORES spectral 
analysis, wind and swell were slightly misaligned. 
PRN08 (or its mirror image) matches very well the 
swell direction and so does PRN 10 along most of 
the track. This result underlines the fact that GNSS- 
R is not sensitive to wind only and that swell has a 
strong impact too. PRN24 has a different behaviour, 
in line with SORES. These two measurements agree 
relatively better with wind direction, although a dis- 
crepancy of 30° is observed at the beginning of the 
track. 

4.3. Slope PDF Isotropy 

It is worth remembering that Elfouhaily's wind- 
driven spectrum predicts a SPI value of 0.65, hardly 
sensitive to wind speed. Here we note that SPI varies 
quite significantly along the track for both GNSS- 
R and SORES. The important departure observed 
from the 0.65 nominal value is probably a signa- 
ture of an under-developed sea and the presence of 
strong swell. Further research should be undertaken 
in order to better understand the potential informa- 
tion contained in this product. 

4.4. Link to wind speed 

On Figure |SJ we have plotted the estimated total 
MSS versus Jason- l's wind speed * together with 

ant wave height, using the algorithm described in [ 10 j. Its stated 
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two models: 

• Elfouhaily's sea-height spectrum, integrated 
for different cut-off wavelengths, and 

• an empirical model proposed by Katzberg 
for L-band, based on a modifica- 
tion of Cox and Munk's relationship: 
MSS=0.9.1(T 3 \/9.48[/ + 6.07C/ 2 , U being 
wind speed (private communication with J.L. 
Garrison, Purdue University). 

We see that both SORES and GNSS-R estimations 
follow Elfouhaily's model trend (MSS obtained by 
integrating the spectrum with the usual cut-off of 3 
times the wavelength) but give higher values of MSS 
(from to 20 to 40% up-shifted). Actually, we have 
found that MSS estimates of PRN08 and 24 are very 
well fitted by Elfouhaily's spectrum with a cut-off of 
one wavelength only. The 20% discrepancy can be 
explained by a strong sea state with a SWH twice as 
high as the one observed during the Cox and Munk's 
experiment (almost 2 m compared to 1 m). At any 
rate, these results indicate that the wind-MSS link 
is not straightforward and that the DMSSa should 
be considered as a self-standing product for oceano- 
graphic users. 

5. Conclusion 

We have reported the first inversion of GNSS- 
R full Delay-Doppler Maps for the retrieval of the 
sea-surface directional mean square slope, DMSSl. 
In addition, we have presented a repetition of the 
Cox and Munk experiment for DMSS opt retrieval 
through inversion of Tilt Azimuth Map of sun glitter 
optical photographs. 

Our results show that both optical and L-band 
total MSS are 20% higher than predicted using 
Elfouhaily's model for the observed wind speed (9 
to 13 m/s). The SPA estimated by GNSS-R matches 
the swell direction with good accuracy for at least 2 
out of 3 PRNs. A new geophysical product has been 
discussed: the slope PDF isotropy, which can be re- 
lated to wind/wave misalignement as well the degree 
of sea development. 

The analysis highlighted the important impact of 
sea-state unmodeled parameters (such as swell) in 
addition to wind stress on the measured DMSSa- 
Since speculometry is sensitive to slope processes 
over a wide range of scales, the link between 
DMSSa and wind is not straightforward: total MSS 



and SPA are definitely affected by swell. Quantita- 
tively, the 20 % bias observed in SORES results can 
be explained by the impact of swell on the elevation 
spectrum. Consequently, DMSSa can and should be 
studied as an independent parameter, of independent 
geophysical value. We note however that the use of 
several wavelenghts could in principle be used to in- 
vert for all the parameters modulating the elevation 
spectrum, a line of future work. 

Let us finally emphasize that the flight was not 
optimized for speculometry (1000 m altitude, 50 m/s 
speed) and that higher/faster flights are needed in the 
future in order to consolidate the concept of DDM 
inversion for DMSSa estimation. 
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Figure 6: Examples of data products and their best-fit models. First column: GNSS-R Delay-Doppler Map (PRN08, arc 01). Second 
column: SORES Tilt-Azimuth Map (photograph 41). First row : Data. Second row : Best-fit Model. Third row : Data-Model 
Residual. 
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Figure 7: DMSS^ estimated along the aircraft track. First row: Total MSS (in dB). Second row: Slope PDF Azimuth. Third row: 
Slope PDF Isotropy. 
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Figure 8: Total MSS versus Jason-l's wind speed. 



